f = oda ~ oda_lag + p_polity2 + gdppc + pop + factor(agency)*factor(year) - 1
f1 = update(f, '.~.+us_troops_2005*independence')
f2 = update(f, '.~.+distance*independence')
f3 = update(f, '.~.+fdi*independence')
f4 = update(f, '.~.+btt*independence')

mod1 = lm(f1, dat)
mod2 = lm(f2, dat)
mod3 = lm(f3, dat)
mod4 = lm(f4, dat)

agency_trends = list(mod1, mod2, mod3, mod4)
se = lapply(agency_trends, function(i) sqrt(diag(sandwich(i))))

print(screenreg(agency_trends, omit.coef='factor', digits=3, override.se=se))
